function y=gmvn(MU,SIGMA,indexConstantValues)

MUreduced=MU(~indexConstantValues);
SIGMAreduced=SIGMA(~indexConstantValues,~indexConstantValues);

randomRealizations=MUreduced + chol(SIGMAreduced)*randn(size(MUreduced));

y=nan(length(MU),1);
y(indexConstantValues)=MU(indexConstantValues);
y(~indexConstantValues)=randomRealizations;

end